Complex Network - Exercises 4¶
Student: Jorge Augusto Salgado Salhani. #USP: 8927418
Student: Nicoly Soares de Almeida. #USP: 8059063
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import random
import pandas as pd
import seaborn as sns
import operator
from tqdm import tqdm
import warnings
warnings.filterwarnings("ignore")
Question 1¶
Here we will analyse the resilience of the network generated models Erdos-Renyi (ER), Barabasi-Albert (BA) and Watts-Strogatz (WS). To accomplish it, we will submit the networks into failures and attacks.
Let us begin with failures. Failures consist on randomly remove a node from the network, and the resilience is measured as how much did the size of the largest component has changed. Firstly, we will take a look at how it does change qualitatively
def failures(H, threshold, print_4instants=False):
"""
FUNCTION:
Submits the graph H into failure process by randomly removing its nodes
PARAMETERS:
H : networkx graph object
threshold: float value that sets the minimum size of the largest component
until the failure process stops
RETURN:
tuple (S, vn): S - vector with the size of the largest component
vn - vector with the amount of nodes that has been taken down
(normalized by the size of the original graph)
"""
fc = fcritical(H)
delta4fc = fc/4
p4i = 0
G = H.copy()
N0 = len(G)
vn = []
S = []
n = 0
i = 0
if print_4instants == True:
while (len(G.nodes()) > int(threshold*N0)):
nodes = np.array(G.nodes)
G.remove_node(random.choice(nodes))
Gcc = sorted(nx.connected_component_subgraphs(G), key=len, reverse=True)
Glc = Gcc[0]
S.append(len(Glc)/N0)
n+=1
vn.append(n/N0)
plt.figure(1, figsize=(15,15))
plt.suptitle('Instant failures before breakdown', fontsize=20)
if n/N0 >= p4i and np.around(p4i, 4) != np.around(fc, 4):
plt.subplot(2,2,i+1)
pos = nx.spring_layout(Glc)
nx.draw_networkx(Glc, node_size=20, with_labels=False, pos=pos, node_color='b')
plt.title('Fraction removed: {0:5.4f}\nLargest Component: {1:5.4f}'.format(p4i, len(Glc)/N0), fontsize=15)
limits = plt.axis('off')
p4i+=delta4fc
i+=1
print(i)
print('Done!')
plt.show(True)
else:
while (len(G.nodes()) > int(threshold*N0)):
nodes = np.array(G.nodes)
G.remove_node(random.choice(nodes))
Gcc = sorted(nx.connected_component_subgraphs(G), key=len, reverse=True)
Glc = Gcc[0]
S.append(len(Glc)/N0)
n+=1
vn.append(n/N0)
return S, vn
Analitically, we can say that exists a fraction of nodes that, after being removed, causes the collapse of the network andthe largest component approaches zero. This measurement is associated to the complexity of a network, as mentioned on previous projects.
The complexity metric is calculated as
$$ \begin{equation} \alpha = \frac{\langle{k^2}\rangle}{\langle{k}\rangle} \end{equation} $$where k is the degree distribution. Therefore, the complexity is the ratio between second and first moment of a distribution. It is related to the robustness/fragility of a network, once it is associated with deviation from mean value. With this quantity, we can also define $f_{c}$ as
$$ \begin{equation} f_c = 1 - \frac{1}{\alpha - 1} \end{equation} $$standing for the fraction of nodes needed to be removed from the network to break it appart. With this result, we can build a function that calculates this quantity
def moment_degree(G, m):
"""
FUNCTION:
Calculates the m-th moment from the degree distribution associated to some given graph G
PARAMETERS:
G: networkx graph object
m: moment degree
RETURN:
float M: m-th moment
"""
M = 0
N = len(G)
for i in G.nodes():
M += G.degree(i)**m
M = M/N
return M
def fcritical(G):
"""
FUNCTION:
Calculates the critical fraction of nodes to breakdown the network after failures occur
PARAMETERS:
G: networkx graph object
RETURN:
f (float): fraction of removed nodes to breakdown
"""
complx = moment_degree(G, 2)
mean = moment_degree(G, 1)
f = 1 - 1/(complx/mean - 1)
return f
Defined our function, we can proceed the failures over the desired networks.
N = 1000
p = 10/N
G = nx.fast_gnp_random_graph(N, p, seed=None, directed=False)
plt.figure(figsize=(15,10))
pos = nx.spring_layout(G)
plt.title('Model Erdos-Renyi', fontsize=20)
nx.draw(G, with_labels = False, node_color='r',
node_size=20, font_size=16, width=.2,pos = pos)
plt.show()
fc = fcritical(G)
Sf, vnf = failures(G, 0.01, print_4instants=True)
N = 1000
av_degree = 10
k = int(av_degree/2)
p = 0.1
G = nx.watts_strogatz_graph(N, k, p, seed=None)
plt.figure(figsize=(15,10))
pos = nx.spring_layout(G)
plt.title('Model Watts-Strogatz with p = 0.1', fontsize=20)
nx.draw(G, with_labels = False, node_color='r',
node_size=20, font_size=16, width=.2,pos = pos)
plt.show()
fc = fcritical(G)
Sf, vnf = failures(G, 0.01, print_4instants=True)
N = 1000
av_degree = 10
m = int(av_degree/2)
G = nx.barabasi_albert_graph(N, m)
plt.figure(figsize=(15,10))
pos = nx.spring_layout(G)
plt.title('Model Barabasi-Albert', fontsize=20)
nx.draw(G, with_labels = False, node_color='r',
node_size=20, font_size=16, width=.2,pos = pos)
plt.show()
fc = fcritical(G)
Sf, vnf = failures(G, 0.01, print_4instants=True)
Now that we saw (visually represented graphs) how the failures alter continuously each network, we can quantify the relation between the largest component and the fraction of removed nodes
def mean_failures_over_n(G, n):
fc = 0
Sf, vnf = failures(G, 0.01)
Sf_tot = np.zeros(len(Sf))
vnf_tot = np.zeros(len(vnf))
for i in range(n):
fc += fcritical(G)/n
Sf, vnf = failures(G, 0.01)
for i in range(len(Sf_tot)):
Sf_tot[i] += Sf[i]/n
vnf_tot[i] += vnf[i]/n
return Sf_tot, vnf_tot, fc
N = 1000
p = 10/N
G = nx.fast_gnp_random_graph(N, p, seed=None, directed=False)
n = 10
Sf_tot_ER, vnf_tot_ER, fc_ER = mean_failures_over_n(G,n)
print('ER Done!')
N = 1000
av_degree = 10
k = int(av_degree/2)
p = 0.1
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Sf_tot_WS, vnf_tot_WS, fc_WS = mean_failures_over_n(G,n)
print('WS Done!')
N = 1000
av_degree = 10
m = int(av_degree/2)
G = nx.barabasi_albert_graph(N, m)
Sf_tot_BA, vnf_tot_BA, fc_BA = mean_failures_over_n(G,n)
print('BA Done!')
plt.figure(figsize=(10,8))
plt.plot(vnf_tot_ER, Sf_tot_ER, '-ob', label='Erdos-Renyi (ER)')
plt.plot(vnf_tot_WS, Sf_tot_WS, '-or', label='Watts-Strogatz (WS)')
plt.plot(vnf_tot_BA, Sf_tot_BA, '-og', label='Barabasi-Albert (BA)')
plt.xlabel('Node fraction removed', fontsize=15)
plt.ylabel('Normalized largest component', fontsize=15)
plt.axvline(x=fc_ER, color='b', label='f critical (ER)')
plt.axvline(x=fc_WS, color='r', label='f critical (WS)')
plt.axvline(x=fc_BA, color='g', label='f critical (BA)')
plt.legend()
plt.grid(True)
plt.show()
Overt $n=10$ samples, it is clear that WS model is the least resilient through failures while ER and BA have near the same behavior. For WS, it will be clearer after we represent how does the resilience change with the rewiring probability $p$
n = 10
N = 1000
av_degree = 10
k = int(av_degree/2)
p = 0.1
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Sf_tot_ER, vnf_tot_ER, fc_ER = mean_failures_over_n(G,n)
print('WS p=0.1 Done!')
p = 0.01
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Sf_tot_WS, vnf_tot_WS, fc_WS = mean_failures_over_n(G,n)
print('WS p=0.01 Done!')
p = 0.001
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Sf_tot_BA, vnf_tot_BA, fc_BA = mean_failures_over_n(G,n)
print('WS p=0.001 Done!')
plt.figure(figsize=(10,8))
plt.plot(vnf_tot_ER, Sf_tot_ER, '-ob', label='WS with p = 0.100')
plt.plot(vnf_tot_WS, Sf_tot_WS, '-or', label='WS with p = 0.010')
plt.plot(vnf_tot_BA, Sf_tot_BA, '-og', label='WS with p = 0.001')
plt.xlabel('Node fraction removed', fontsize=15)
plt.ylabel('Normalized largest component', fontsize=15)
plt.axvline(x=fc_ER, color='b', label='f critical (WS, p=0.100)')
plt.axvline(x=fc_WS, color='r', label='f critical (WS, p=0.010)')
plt.axvline(x=fc_BA, color='g', label='f critical (WS, p=0.001')
plt.legend()
plt.grid(True)
plt.show()
As we can see, also taking the average over $n=10$ samples, the higher the value of probability $p$, the closer the shape resembles the ER model. It happens because for $p=1$ the rewiring happens for each edge, thus connecting the whole network as a random generated one (closer to ER). On the other extreme, for $p=0$ the network is linearly shaped and the removal of a few nodes suffices to break it appart.
Now we can perform the same analysis for attacks. Attacks consist on remove a node from the network based on its degree, meaning that the higher the degree, the higher the node's degree, the quicker it will be removed from the network.
Forthis purpouse, we can start building a function that returns the most connected node and, after, its exclusion
def most_connected_node(G):
"""
FUNCTION:
Get the node with the maximum degree over the network
PARAMETERS:
G: networkx graph object
RETURNS:
node (int): number of the node with the highest degree
"""
maxk = 0
node = 0
for i in G.nodes():
if G.degree(i) >= maxk:
maxk = G.degree(i)
node = i
return node
def attacks(H, threshold, print_4instants=False, fc_by_hand=0):
"""
FUNCTION:
Submits the graph H into atack process by removing the node with highest degree
PARAMETERS:
H : networkx graph object
threshold: float value that sets the minimum size of the largest component
until the failure process stops
RETURN:
tuple (S, vn): S - vector with the size of the largest component
vn - vector with the amount of nodes that has been taken down
(normalized by the size of the original graph)
"""
G = H.copy()
N0 = len(G)
vn = []
S = []
n = 0
i = 0
if print_4instants == True:
if fc_by_hand==0:
fc = fcritical(G)
else:
fc = fc_by_hand
delta4fc = fc/4
p4i = 0
while (len(G.nodes()) > int(threshold*N0)):
node = most_connected_node(G)
G.remove_node(node)
Gcc = sorted(nx.connected_component_subgraphs(G), key=len, reverse=True)
Glc = Gcc[0]
S.append(len(Glc)/N0)
n+=1
vn.append(n/N0)
plt.figure(1, figsize=(15,15))
plt.suptitle('Instant attacks before breakdown', fontsize=20)
if n/N0 >= p4i and np.around(p4i, 4) != np.around(fc, 4):
plt.subplot(2,2,i+1)
pos = nx.spring_layout(Glc)
nx.draw_networkx(Glc, node_size=20, with_labels=False, pos=pos, node_color='b')
plt.title('Fraction removed: {0:5.4f}\nLargest Component: {1:5.4f}'.format(p4i, len(Glc)/N0), fontsize=15)
limits = plt.axis('off')
p4i+=delta4fc
i+=1
print(i)
print('Done!')
plt.show(True)
else:
while (len(G.nodes()) > int(threshold*N0)):
node = most_connected_node(G)
G.remove_node(node)
Gcc = sorted(nx.connected_component_subgraphs(G), key=len, reverse=True)
Glc = Gcc[0]
S.append(len(Glc)/N0)
n+=1
vn.append(n/N0)
return S, vn
def mean_attacks_over_n(G, n):
Sf, vnf = attacks(G, 0.01)
Sf_tot = np.zeros(len(Sf))
vnf_tot = np.zeros(len(vnf))
for i in range(n):
Sf, vnf = attacks(G, 0.01)
for i in range(len(Sf_tot)):
Sf_tot[i] += Sf[i]/n
vnf_tot[i] += vnf[i]/n
return Sf_tot, vnf_tot
Defined our functions, we submit the networks under attacks
N = 1000
p = 10/N
G = nx.fast_gnp_random_graph(N, p, seed=None, directed=False)
plt.figure(figsize=(15,10))
pos = nx.spring_layout(G)
plt.title('Model Erdos-Renyi', fontsize=20)
nx.draw(G, with_labels = False, node_color='r',
node_size=20, font_size=16, width=.2,pos = pos)
plt.show()
Sf, vnf = attacks(G, 0.01, print_4instants=True, fc_by_hand=0.6)
N = 1000
av_degree = 10
k = int(av_degree/2)
p = 0.1
G = nx.watts_strogatz_graph(N, k, p, seed=None)
plt.figure(figsize=(15,10))
pos = nx.spring_layout(G)
plt.title('Model Watts-Strogatz with p = 0.1', fontsize=20)
nx.draw(G, with_labels = False, node_color='r',
node_size=20, font_size=16, width=.2,pos = pos)
plt.show()
Sf, vnf = attacks(G, 0.01, print_4instants=True, fc_by_hand=0.3)
N = 1000
av_degree = 10
m = int(av_degree/2)
G = nx.barabasi_albert_graph(N, m)
plt.figure(figsize=(15,10))
pos = nx.spring_layout(G)
plt.title('Model Barabasi-Albert', fontsize=20)
nx.draw(G, with_labels = False, node_color='r',
node_size=20, font_size=16, width=.2,pos = pos)
plt.show()
Sf, vnf = attacks(G, 0.01, print_4instants=True,fc_by_hand=0.4)
Similarly as before, we can make a quantitative analysis for resilience of each network
N = 1000
p = 10/N
G = nx.fast_gnp_random_graph(N, p, seed=None, directed=False)
n = 10
Sf_tot_ER, vnf_tot_ER = mean_attacks_over_n(G,n)
print('ER Done!')
N = 1000
av_degree = 10
k = int(av_degree/2)
p = 0.1
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Sf_tot_WS, vnf_tot_WS = mean_attacks_over_n(G,n)
print('WS Done!')
N = 1000
av_degree = 10
m = int(av_degree/2)
G = nx.barabasi_albert_graph(N, m)
Sf_tot_BA, vnf_tot_BA = mean_attacks_over_n(G,n)
print('BA Done!')
plt.figure(figsize=(10,8))
plt.plot(vnf_tot_ER, Sf_tot_ER, '-ob', label='Erdos-Renyi (ER)')
plt.plot(vnf_tot_WS, Sf_tot_WS, '-or', label='Watts-Strogatz (WS)')
plt.plot(vnf_tot_BA, Sf_tot_BA, '-og', label='Barabasi-Albert (BA)')
plt.xlabel('Node fraction removed', fontsize=15)
plt.ylabel('Normalized largest component', fontsize=15)
plt.legend()
plt.grid(True)
plt.show()
For attacks, distinctly from failures, ER outperforms the BA model. The WS will be analyse further, with varying the rewiring probability $p$. Once BA model presents some highly connected nodes, the attacks are more effective her than on ER, tearing the network appart easily by removing the biggest hubs at the initial times
Finally, for different probabilitis for Watts-Strogatz network model
n = 10
N = 1000
av_degree = 10
k = int(av_degree/2)
p = 0.1
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Sf_tot_ER, vnf_tot_ER = mean_attacks_over_n(G,n)
print('WS p=0.1 Done!')
p = 0.01
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Sf_tot_WS, vnf_tot_WS = mean_attacks_over_n(G,n)
print('WS p=0.01 Done!')
p = 0.001
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Sf_tot_BA, vnf_tot_BA = mean_attacks_over_n(G,n)
print('WS p=0.001 Done!')
plt.figure(figsize=(10,8))
plt.plot(vnf_tot_ER, Sf_tot_ER, '-ob', label='WS with p = 0.100')
plt.plot(vnf_tot_WS, Sf_tot_WS, '-or', label='WS with p = 0.010')
plt.plot(vnf_tot_BA, Sf_tot_BA, '-og', label='WS with p = 0.001')
plt.xlabel('Node fraction removed', fontsize=15)
plt.ylabel('Normalized largest component', fontsize=15)
plt.legend()
plt.grid(True)
plt.show()
Interestingly, the altering on $p$ value starts closer to the behavior of WS submitted to failures, for $p = 0.001$, it goes throughout a highly susceptible to breakdown as $p$ increases to $p = 0.01$ and, again, returns to a resilient mode for $p=0.1$, which is closer to the behavior with $p = 0.001$. Therefore we can conclude that there exists a value for $p$ that minimizes the network resilience
def find(v, i):
"""
FUNCTION:
Found the position of a element i in a given vector v
PARAMETERS:
v: vector to be analysed
i: element to be found in vector v
RETURN:
vector: Containing the position of the element i in vector v (if i is in v)
"""
l = []
pos = 0
for x in v:
if (x == i):
l.append(pos)
pos+=1
return l
P = np.arange(0.001, 0.1, 0.005)
n = 10
N = 1000
av_degree = 10
k = int(av_degree/2)
min_grad = np.zeros(len(P))
T = 10
for i in range(T):
count = 0
for p in P:
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Sf_tot_ER, vnf_tot_ER = mean_attacks_over_n(G,n)
grad = np.gradient(Sf_tot_ER)
pos = find(grad, min(grad))
min_grad[count] += vnf_tot_ER[pos[0]]/T
count+=1
plt.figure(figsize=(8,6))
plt.title('Highest Sudden Change Before Breakdown for WS Through Attacks ', fontsize=15)
plt.xlabel('Rewiring probability p', fontsize=13)
plt.ylabel('Fraction of removed nodes (abrupt changes)', fontsize=13)
plt.plot(P, min_grad, '-o')
plt.grid(True)
plt.show()
On this last graph we analyse how does the fraction of removed nodes at the inflection point on graph measuring the 'largest component size' $\textit{versus}$ 'fraction of removed nodes' varies as the rewiring probability $p$ changes. As previously mentioned, the resilience on WS network does not linearly correlates with the value of $p$, but presents transitory behavior where the resilience starts high with $p \approx 0$, decreases and starts again to increase. Therefore on this last graph we tried to quantify how does the value of $p$ correlates to the resilience of WS model. It turns out that $p \approx 0.3$ provides the lowest resilience throught attacks.
Considering now some real networks, we can analyse their resilience in terms of attacks and failures, so we must have some functions to read networkx graphs
def read_gml_and_Gcc(path, nodetype=int):
"""
FUNCTION:
Read a file gml type that contains node connections
and their weights as float.
PARAMETERS:
path: file path as string located in computer
RETURN:
tuple: (pos, G), representing the position object
(networkx as nx, spring_layout()) for draw purpose and
the Graph G itself made from file
"""
G = nx.read_gml(path)
G = nx.convert_node_labels_to_integers(G, first_label=0)
labels = G.nodes()
pos = nx.spring_layout(G)
G = G.to_undirected()
G.remove_edges_from(G.selfloop_edges())
Gcc = sorted(nx.connected_component_subgraphs(G), key=len, reverse=True)
G = Gcc[0]
return pos, G
def read_txt_and_Gcc(path, comment='%', nodetype=int):
"""
FUNCTION:
Read a file that contains node connections and their
weights as float. Example of readable file
---------
%This is an example file.
%NodeA NodeB Weight
0 1 0.5
1 3 2.0
2 3 1.2
---------
Comments can be marked as wish ('%' is set as default)
PARAMETERS:
path: file path as string located in computer
comment: char (or string) used in file to mark comments
RETURN:
tuple: (pos, G), representing the position object
(networkx as nx, spring_layout()) for draw purpose and
the Graph G itself made from file
"""
G = nx.read_edgelist(path, comments=comment, nodetype=nodetype, data=(('weight', float),))
G = nx.convert_node_labels_to_integers(G, first_label=0)
labels = G.nodes()
pos = nx.spring_layout(G)
G = G.to_undirected()
G.remove_edges_from(G.selfloop_edges())
# Gcc = sorted(nx.connected_component_subgraphs(G), key=len, reverse=True)
# G = Gcc[0]
return pos, G
def read_csv_and_Gcc(path, source='u', target='v', nodetype=int):
"""
FUNCTION:
Read a file gml type that contains node connections
and their weights as float.
PARAMETERS:
path: file path as string located in computer
RETURN:
tuple: (pos, G), representing the position object
(networkx as nx, spring_layout()) for draw purpose and
the Graph G itself made from file
"""
import pandas as pd
df = pd.read_csv(path)
df.head()
G = nx.from_pandas_edgelist(df, source=source, target=target)
G = nx.convert_node_labels_to_integers(G, first_label=0)
labels = G.nodes()
pos = nx.spring_layout(G)
G = G.to_undirected()
G.remove_edges_from(G.selfloop_edges())
Gcc = sorted(nx.connected_component_subgraphs(G), key=len, reverse=True)
G = Gcc[0]
return pos, G
pos_protein, G_protein = read_txt_and_Gcc('/home/jorge/Documents/redes_complexas/data/ppi.maayan-vidal.txt')
pos_celegans7, G_celegans7 = read_txt_and_Gcc('/home/jorge/Documents/redes_complexas/wi2007.txt', comment='#', nodetype=str)
pos_mosquito, G_mosquito = read_csv_and_Gcc('/home/jorge/Documents/redes_complexas/pntd.0005879.s001.csv', source='Dengue Protein', target='Human Protein', nodetype=str)
For visualisation sake, we can show those networks
plt.figure(figsize=(15,10))
plt.title('Human Proteome', fontsize=20)
nx.draw(G_protein, with_labels = False, node_color='r',
node_size=15, font_size=16, width=.2, pos=pos_protein)
plt.show()
plt.figure(figsize=(15,10))
plt.title('C. elegans neural network', fontsize=20)
nx.draw(G_celegans7, with_labels = False, node_color='r',
node_size=15, font_size=16, width=.2, pos=pos_celegans7)
plt.show()
plt.figure(figsize=(15,10))
plt.title('Human-Dengue Proteins Interactome', fontsize=20)
nx.draw(G_mosquito, with_labels = False, node_color='r',
node_size=30, font_size=16, width=.2, pos=pos_mosquito)
plt.show()
Now we can submit those networks into failures and attacks. As before, we can play with the visual breaks of the networks before breakdown and the quantitativity analysis
fc = fcritical(G_protein)
Sf, vnf = failures(G_protein, 0.01, print_4instants=True)
fc = fcritical(G_celegans7)
Sf, vnf = failures(G_celegans7, 0.01, print_4instants=True)
fc = fcritical(G_mosquito)
Sf, vnf = failures(G_mosquito, 0.01, print_4instants=True)
Now, under attacks
Sf, vnf = attacks(G_protein, 0.01, print_4instants=True, fc_by_hand=0.2)
Sf, vnf = attacks(G_celegans7, 0.01, print_4instants=True, fc_by_hand=0.1)
Sf, vnf = attacks(G_mosquito, 0.01, print_4instants=True, fc_by_hand=0.03)
Now, for quantitative analysis
plt.figure(1, figsize=(16,6))
plt.suptitle('Resilience of Real Networks', fontsize=20)
plt.subplot(1,3,1)
plt.title('Human Proteome', fontsize=15)
fc = fcritical(G_protein)
Sf, vnf = failures(G_protein, 0.01)
Sa, vna = attacks(G_protein, 0.01)
plt.plot(vnf, Sf, '-or', label='Failures')
plt.plot(vna, Sa, '-ob', label='Attacks')
plt.ylabel('Size largest component', fontsize=13)
plt.axvline(x=fc, color='g', label='fc failures')
plt.legend()
plt.subplot(1,3,2)
plt.title('C. elegans Neurals', fontsize=15)
fc = fcritical(G_celegans7)
Sf, vnf = failures(G_celegans7, 0.01)
Sa, vna = attacks(G_celegans7, 0.01)
plt.plot(vnf, Sf, '-or', label='Failures')
plt.plot(vna, Sa, '-ob', label='Attacks')
plt.xlabel('Fraction removed nodes', fontsize=13)
plt.axvline(x=fc, color='g', label='fc failures')
plt.legend()
plt.subplot(1,3,3)
plt.title('Human-Dengue Interactome', fontsize=15)
fc = fcritical(G_mosquito)
Sf, vnf = failures(G_mosquito, 0.01)
Sa, vna = attacks(G_mosquito, 0.01)
plt.plot(vnf, Sf, '-or', label='Failures')
plt.plot(vna, Sa, '-ob', label='Attacks')
plt.axvline(x=fc, color='g', label='fc failures')
plt.legend()
plt.show()
Comparing the real networks presented above, we can conclude that the Human proteome happens to have the most and the C. elegans, the least, resilience throughout failures. On the other hand, throughout attacks thehuman proteome is also the most resilient and the human-dengue interactome, the least.
Visually we could conclude that human-dengue interactome should be the least resilient, once it has the most apparent hubs built from a single node, however the human interactome resilience could be dues to its highest amount of nodes when compared to the others, which provides resilience to the network besides its topology
Considering the Fortunato benchmark for community organization, we can also see how it is resilient through failures and attacks
from networkx.algorithms.community import LFR_benchmark_graph
N = 128
tau1 = 3
tau2 = 1.5
k =16
minc = 32
maxc = 32
plt.figure(1, figsize=(15,6))
plt.suptitle('Fortunato benchmark networks', fontsize=20)
plt.subplot(1,3,1)
mu = 0.1
plt.title('Mixing $\mu$ = 0.1', fontsize=15)
G0 = LFR_benchmark_graph(n = N, tau1 = tau1, tau2 = tau2, mu = mu, min_degree = k,
max_degree = k, min_community=minc, max_community = maxc, seed = 10)
pos=nx.spring_layout(G0)
nx.draw(G0, with_labels = False, nodecolor='r', node_size=50, font_size=16, width=1, pos=pos)
plt.subplot(1,3,2)
plt.title('Mixing $\mu$ = 0.3', fontsize=15)
mu = 0.3
G1 = LFR_benchmark_graph(n = N, tau1 = tau1, tau2 = tau2, mu = mu, min_degree = k,
max_degree = k, min_community=minc, max_community = maxc, seed = 10)
pos=nx.spring_layout(G1)
nx.draw(G1, with_labels = False, nodecolor='r', node_size=50, font_size=16, width=1, pos=pos)
plt.subplot(1,3,3)
plt.title('Mixing $\mu$ = 0.5', fontsize=15)
mu = 0.5
G2 = LFR_benchmark_graph(n = N, tau1 = tau1, tau2 = tau2, mu = mu, min_degree = k,
max_degree = k, min_community=minc, max_community = maxc, seed = 10)
pos=nx.spring_layout(G2)
nx.draw(G2, with_labels = False, nodecolor='r', node_size=50, font_size=16, width=1, pos=pos)
plt.show(True)
fc = fcritical(G0)
Sf, vnf = failures(G0, 0.01, print_4instants=True)
fc = fcritical(G1)
Sf, vnf = failures(G1, 0.01, print_4instants=True)
fc = fcritical(G2)
Sf, vnf = failures(G2, 0.01, print_4instants=True)
fc = fcritical(G0)
Sf, vnf = attacks(G0, 0.01, print_4instants=True)
fc = fcritical(G1)
Sf, vnf = attacks(G1, 0.01, print_4instants=True)
fc = fcritical(G2)
Sf, vnf = attacks(G2, 0.01, print_4instants=True)
plt.figure(1, figsize=(16,6))
plt.suptitle('Resilience of Fortunato benchmark', fontsize=20)
plt.subplot(1,3,1)
plt.title('Mixing $\mu$ = 0.1', fontsize=15)
fc = fcritical(G0)
Sf, vnf = failures(G0, 0.01)
Sa, vna = attacks(G0, 0.01)
plt.plot(vnf, Sf, '-or', label='Failures')
plt.plot(vna, Sa, '-ob', label='Attacks')
plt.ylabel('Size largest component', fontsize=13)
plt.axvline(x=fc, color='g', label='fc failures')
plt.legend()
plt.subplot(1,3,2)
plt.title('Mixing $\mu$ = 0.3', fontsize=15)
fc = fcritical(G1)
Sf, vnf = failures(G1, 0.01)
Sa, vna = attacks(G1, 0.01)
plt.plot(vnf, Sf, '-or', label='Failures')
plt.plot(vna, Sa, '-ob', label='Attacks')
plt.xlabel('Fraction removed nodes', fontsize=13)
plt.axvline(x=fc, color='g', label='fc failures')
plt.legend()
plt.subplot(1,3,3)
plt.title('Mixing $\mu$ = 0.5', fontsize=15)
fc = fcritical(G2)
Sf, vnf = failures(G2, 0.01)
Sa, vna = attacks(G2, 0.01)
plt.plot(vnf, Sf, '-or', label='Failures')
plt.plot(vna, Sa, '-ob', label='Attacks')
plt.axvline(x=fc, color='g', label='fc failures')
plt.legend()
plt.show()
The Fortunato benchmark happens to have no huge distinction on its resilience when varying the mixing parameter. The failures and attacks seems to breakdown the network at approximatelly the same fraction of nodes being removed. The only distinction between them relates to the smoothness of the transition between connected and broken. While increasing the value of mixing $\mu$, the network transits more gradually between those two states, which does make sense, since the highest the value of $\mu$, the closer the network resembles a random connected one, which dilutes the hubs and, therefore, softens the transition from a connected to a broken network
One of the algorithms that helps the study of assortativenness of networks is the Xulvi-Brunet-Sokolov algorithm. It increases or decreases the assortativenness while keep the network degree distribution characteristic unchanged. The algorithm can be described as follows:
Choose two random edges from the graph ($e_1$, $e_2$)
Check the degree of each one of the 4 nodes and sort them ($k_1$, $k_2$, $k_3$, $k_4$, with $k_1 > k_2 > k_3 > k_4$)
Increase the assortativenness: Connects $k_1$ and $k_2$, $k_3$ and $k_4$, breaking the previous connection between them
Decrease the assortativenness: Connects $k_1$ and $k_4$, $k_2$ and $k_3$, breaking the previous connection between them
Below we build a function for each of the cases described at 3. and 4.
def xulvi_brunet_sokolov_assortative(G, p):
edges = list(G.edges())
chosen1 = np.random.randint(0, len(edges)-1)
e1 = edges[chosen1]
chosen2 = np.random.randint(0, len(edges)-1)
e2 = edges[chosen2]
k1 = len(list(G.neighbors(e1[0])))
k2 = len(list(G.neighbors(e1[1])))
k3 = len(list(G.neighbors(e2[0])))
k4 = len(list(G.neighbors(e2[1])))
#print(e1, e2)
ks = {e1[0]:k1, e1[1]:k2, e2[0]:k3, e2[1]:k4}
p = 1
if len(ks) == 4:
#print(ks)
ks_sort = sorted(ks.items(), key=operator.itemgetter(1), reverse=True)
#print(ks_sort)
#print(G.has_edge(ks_sort[0][0], ks_sort[1][0]), G.has_edge(ks_sort[2][0], ks_sort[3][0]))
if G.has_edge(ks_sort[0][0], ks_sort[1][0]) == False and np.random.rand() < p:
G.remove_edge(e1[0],e1[1])
G.add_edge(ks_sort[0][0], ks_sort[1][0])
G.remove_edge(e2[0],e2[1])
G.add_edge(ks_sort[2][0], ks_sort[3][0])
#print(G.has_edge(ks_sort[0][0], ks_sort[1][0]), G.has_edge(ks_sort[2][0], ks_sort[3][0]))
def xulvi_brunet_sokolov_disassortative(G, p):
edges = list(G.edges())
chosen1 = np.random.randint(0, len(edges)-1)
e1 = edges[chosen1]
chosen2 = np.random.randint(0, len(edges)-1)
e2 = edges[chosen2]
k1 = len(list(G.neighbors(e1[0])))
k2 = len(list(G.neighbors(e1[1])))
k3 = len(list(G.neighbors(e2[0])))
k4 = len(list(G.neighbors(e2[1])))
#print(e1, e2)
ks = {e1[0]:k1, e1[1]:k2, e2[0]:k3, e2[1]:k4}
p = 1
if len(ks) == 4:
#print(ks)
ks_sort = sorted(ks.items(), key=operator.itemgetter(1), reverse=True)
#print(ks_sort)
#print(G.has_edge(ks_sort[0][0], ks_sort[1][0]), G.has_edge(ks_sort[2][0], ks_sort[3][0]))
if G.has_edge(ks_sort[0][0], ks_sort[3][0]) == False and np.random.rand() < p:
G.remove_edge(e1[0],e1[1])
G.add_edge(ks_sort[0][0], ks_sort[3][0])
G.remove_edge(e2[0],e2[1])
G.add_edge(ks_sort[1][0], ks_sort[2][0])
#print(G.has_edge(ks_sort[0][0], ks_sort[1][0]), G.has_edge(ks_sort[2][0], ks_sort[3][0]))
Now, provided a BA network, we can see how does the algorithm changes the assortativenness of the network
N = 500
av_degree = 10
m = int(av_degree/2)
G = nx.barabasi_albert_graph(N, m)
p = 1
T = 15000
plt.figure(figsize=(15,12))
plt.suptitle('Xulvi-Brunet-Sokolov Assortative', fontsize=20)
count = 1
for i in range(T+1):
xulvi_brunet_sokolov_assortative(G, p)
if i%5000 == 0:
#print(count)
plt.subplot(2,2,count)
plt.axis('off')
pos = nx.spring_layout(G)
nx.draw_networkx(G, node_size=50, pos=pos, with_labels=False)
assort = nx.degree_assortativity_coefficient(G)
plt.title('Assortativity: {0:.4f}'.format(assort), fontsize=15)
count+=1
#print('Done!')
plt.show()
N = 500
av_degree = 10
m = int(av_degree/2)
G = nx.barabasi_albert_graph(N, m)
p = 1
T = 15000
plt.figure(figsize=(15,12))
plt.suptitle('Xulvi-Brunet-Sokolov Disassortative', fontsize=20)
count = 1
for i in range(T+1):
xulvi_brunet_sokolov_disassortative(G, p)
if i%5000 == 0:
#print(count)
plt.subplot(2,2,count)
plt.axis('off')
pos = nx.spring_layout(G)
nx.draw_networkx(G, node_size=50, pos=pos, with_labels=False)
assort = nx.degree_assortativity_coefficient(G)
plt.title('Assortativity: {0:.4f}'.format(assort), fontsize=13)
count+=1
#print('Done!')
plt.show()
Once stabilished the possibility of change on assortativenness, we can generate a BA network, reduce its assortativennes and start to measure how its resilience change while the assortative increases
plt.figure(1, figsize=(15,12))
plt.suptitle('Resilience for Barabasi-Albert with Assortativity', fontsize=20)
N = 500
av_degree = 10
m = int(av_degree/2)
G = nx.barabasi_albert_graph(N, m)
p = 1
T = 15000
count = 1
for i in range(T+1):
xulvi_brunet_sokolov_disassortative(G, p)
T = 10000
for i in range(T+1):
xulvi_brunet_sokolov_assortative(G, p)
assort = nx.degree_assortativity_coefficient(G)
if i == 0 or i == 500 or i == 5000 or i ==10000:
Sf, vnf = failures(G, 0.01)
Sa, vna = attacks(G, 0.01)
#print(count)
plt.subplot(2,2,count)
count+=1
fc = fcritical(G)
plt.title('Assortativity: {0:.4f}'.format(assort), fontsize=15)
plt.plot(vnf, Sf, '-or', label='Failures')
plt.plot(vna, Sa, '-ob', label='Attacks')
plt.axvline(x=fc, color='g', label='fc failures')
plt.legend()
if i == 0 or i == 5000:
plt.ylabel('Size Largest Component', fontsize=13)
if i == 5000 or i == 10000:
plt.xlabel('Fraction removed nodes', fontsize=13)
#print('Done!')
plt.show()
It can be concluded that assortativenness is directly correlated to resilience when networks are submitted to attacks (for failures, there are no expressive changes). For disassortative networks (closer to -1, up-left plot), nodes with higher degrees are linked to the lowest degree ones, therefore we have a network that is similar to the human-dengue interactome, which is very susceptible to breakdown through attacks. Meanwhile for assortative networks (closer to +1, down-right plot) higher-degree nodes also are mostly connected with high-degree nodes, which makes the network less susceptible to breakdown, since the attacks do not automatically reduces the largest component of the graph.
Question 2¶
In this question we will attend the epidemic spreading on different networks for Erdos-Renyi, Watts-Strogatz and Barabasi-Albert models.
There are two main models we can set to evolve on each network. They are called SIS (Susceptible-Infected-Suscetible) and SIR (Susceptible-Infected-Resistent)
The SIS model consists on individuals (represented by each node on a graph) that can be at state $S$ and, when connected to someone who is at state $I$, can also become $I$ with probability $\beta \in [0,1]$. After the next steps, it can return to state $S$ with probability $\mu \in [0,1]$.
It can be displayed as:
$$ \begin{equation} S \xrightarrow{\beta} I \xrightarrow{\mu} S \end{equation} $$The SIR model consists on individuals that can be at state $S$ and, when in contact to someone who is at state $I$, can also become $I$ with probability $\beta$. After the next steps, it can become $R$ (with no chances to be infected again) with probability $\mu$
It can be displayed as:
$$ \begin{equation} S \xrightarrow{\beta} I \xrightarrow{\mu} R \end{equation} $$At continuous time, the dynamics can be written in form of differential equation such as, for SIS model,
$$ \begin{align} \frac{d S}{d t} &= -\beta S I + \mu I \\ \frac{d I}{d t} &= \beta S I - \mu I \end{align} $$and, for SIR
$$ \begin{align} \frac{d S}{d t} &= -\beta S I \\ \frac{d I}{d t} &= \beta S I - \mu I \\ \frac{d R}{d t} &= \mu I \end{align} $$Now we can build some function in order to perform the dynamics over a given network
def find(v, i):
"""
FUNCTION:
Found the position of a element i in a given vector v
PARAMETERS:
v: vector to be analysed
i: element to be found in vector v
RETURN:
vector: Containing the position of the element i in vector v (if i is in v)
"""
l = []
pos = 0
for x in v:
if (x == i):
l.append(pos)
pos+=1
return l
def SIR_model_spread(G, beta, mu, Tmax):
"""
FUNCTION:
Performs the SIR dynamic process on a given graph, placing a infectious agent (seed)
on each node each time and spreading it over the others
PARAMETERS:
G: networkx graph object
beta: probability of infection of a Susceptible agent
mu: probability of recovery of a Infected agent
Tmax: Maximum time over which the simulation will proceed
RETURN:
tuple (vt, av_infect): Contains vectors vt: time span
av_infect: average number of infectious nodes
"""
N = len(G)
av_infec = np.zeros(Tmax)
for seed_node in G.nodes():
vector_states = np.zeros(N)
vector_states[seed_node] = 1
ninfected = 1
t = 0
infected = list()
vt = list()
infect = list()
while ninfected > 0 and t < Tmax:
infected = find(vector_states, 1)
for i in infected:
neigs = G.neighbors(i)
for j in neigs:
if np.random.rand() < beta:
if (vector_states[j] != 2):
vector_states[j] = 1
for k in infected:
if np.random.rand() < mu:
vector_states[k] = 2
ninfected = len(find(vector_states, 1))
infect.append(ninfected/N)
t+=1
vt.append(t)
for x in np.arange(0, len(infect)):
av_infec[x] = av_infec[x] + infect[x]
av_infec = av_infec/len(G.nodes())
vt = np.arange(0, Tmax)
return vt, av_infec
def SIS_model_spread(G, beta, mu, Tmax):
"""
FUNCTION:
Performs the SIS dynamic process on a given graph, placing a infectious agent (seed)
on each node each time and spreading it over the others
PARAMETERS:
G: networkx graph object
beta: probability of infection of a Susceptible agent
mu: probability of recovery of a Infected agent
Tmax: Maximum time over which the simulation will proceed
RETURN:
tuple (vt, av_infect): Contains vectors vt: time span
av_infect: average number of infectious nodes
"""
N = len(G)
av_infec = np.zeros(Tmax)
for seed_node in G.nodes():
vector_states = np.zeros(N)
vector_states[seed_node] = 1
ninfected = 1
t = 0
infected = list()
vt = list()
infect = list()
for t in range(Tmax):
infected = find(vector_states, 1)
for i in infected:
neigs = G.neighbors(i)
for j in neigs:
if np.random.rand() < beta:
vector_states[j] = 1
for k in infected:
if np.random.rand() < mu:
vector_states[k] = 0
ninfected = len(find(vector_states, 1))
infect.append(ninfected/N)
t+=1
vt.append(t)
for x in np.arange(0, len(infect)):
av_infec[x] = av_infec[x] + infect[x]
av_infec = av_infec/len(G.nodes())
vt = np.arange(0, Tmax)
return vt, av_infec
plt.figure(1, figsize=(15,10))
plt.suptitle('SIS and SIR spread models (Infected only)', fontsize=20)
####################
plt.subplot(2,3,1)
plt.title('Erdos-Renyi', fontsize=15)
N = 500
av_degree = 8
p = 10/N
G = nx.fast_gnp_random_graph(N, p, seed=None, directed=False)
beta = 0.3
mu = 1
Tmax = 20
vt, inf = SIR_model_spread(G, beta, mu, Tmax)
plt.plot(vt, inf, 'bo-', label='SIR')
vt, inf = SIS_model_spread(G, beta, mu, Tmax)
plt.plot(vt, inf, 'ro-', label='SIS')
plt.ylabel('Infected nodes', fontsize=13)
plt.legend()
###################
plt.subplot(2,3,3)
plt.title('Barabasi-Albert', fontsize=15)
m = int(av_degree/2)
G = nx.barabasi_albert_graph(N, m)
vt, inf = SIR_model_spread(G, beta, mu, Tmax)
plt.plot(vt, inf, 'bo-', label='SIR')
vt, inf = SIS_model_spread(G, beta, mu, Tmax)
plt.plot(vt, inf, 'ro-', label='SIS')
plt.legend()
###################
plt.subplot(2,3,4)
plt.title('Watts-Strogatz p = 0.1', fontsize=15)
k = int(av_degree/2)
p = 0.1
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Tmax = 40
vt, inf = SIR_model_spread(G, beta, mu, Tmax)
plt.plot(vt, inf, 'bo-', label='SIR')
vt, inf = SIS_model_spread(G, beta, mu, Tmax)
plt.plot(vt, inf, 'ro-', label='SIS')
plt.ylabel('Infected nodes', fontsize=13)
plt.legend()
####################
plt.subplot(2,3,5)
plt.title('Watts-Strogatz p = 0.01', fontsize=15)
k = int(av_degree/2)
p = 0.01
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Tmax = 40
vt, inf = SIR_model_spread(G, beta, mu, Tmax)
plt.plot(vt, inf, 'bo-', label='SIR')
vt, inf = SIS_model_spread(G, beta, mu, Tmax)
plt.xlabel('Time', fontsize=13)
plt.plot(vt, inf, 'ro-', label='SIS')
plt.legend()
####################
plt.subplot(2,3,6)
plt.title('Watts-Strogatz p = 0.001', fontsize=15)
k = int(av_degree/2)
p = 0.001
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Tmax = 40
vt, inf = SIR_model_spread(G, beta, mu, Tmax)
plt.plot(vt, inf, 'bo-', label='SIR')
vt, inf = SIS_model_spread(G, beta, mu, Tmax)
plt.plot(vt, inf, 'ro-', label='SIS')
plt.legend()
plt.show()
Finally we can consider the critical probability, namely $\lambda_c$ such that $\lambda = \beta / \mu$, that is the ratio between the infectious and recovery rate of a disease.
def critical_SIS_SIR_models(G, mu, Tmax):
mu = 1
vlabda = list()
rhoI_SIS = list()
rhoI_SIR = list()
for lbd in np.arange(0, 1, 0.05):
beta = lbd*mu
vt, av_infec = SIS_model_spread(G, beta, mu, Tmax)
vlabda.append(lbd)
rhoI_SIS.append(np.mean(av_infec))
vt, av_infec = SIR_model_spread(G, beta, mu, Tmax)
#print(av_infec)
rhoI_SIR.append(np.mean(av_infec))
#print(lbd)
return vlabda, rhoI_SIS, rhoI_SIR
plt.figure(1, figsize=(15,10))
plt.suptitle('SIS and SIR spread models (Infected only)', fontsize=20)
####################
plt.subplot(2,3,1)
plt.title('Erdos-Renyi', fontsize=15)
N = 500
p = 10/N
G = nx.fast_gnp_random_graph(N, p, seed=None, directed=False)
mu = 1
Tmax = 20
vlabda, rhoI_SIS, rhoI_SIR = critical_SIS_SIR_models(G, mu, Tmax)
plt.plot(vlabda, rhoI_SIS, 'ro--', label='SIS')
plt.plot(vlabda, rhoI_SIR, 'bo--', label='SIR')
plt.ylabel('Infected nodes', fontsize=13)
lbdac = moment_degree(G, 1)/moment_degree(G,2)
plt.axvline(lbdac, color='g', label='$\lambda_c$ = {0:5.3f}'.format(lbdac))
plt.legend()
###################
plt.subplot(2,3,3)
plt.title('Barabasi-Albert', fontsize=15)
m = int(av_degree/2)
G = nx.barabasi_albert_graph(N, m)
vlabda, rhoI_SIS, rhoI_SIR = critical_SIS_SIR_models(G, mu, Tmax)
plt.plot(vlabda, rhoI_SIS, 'ro--', label='SIS')
plt.plot(vlabda, rhoI_SIR, 'bo--', label='SIR')
lbdac = moment_degree(G, 1)/moment_degree(G,2)
plt.axvline(lbdac, color='g', label='$\lambda_c$ = {0:5.3f}'.format(lbdac))
plt.legend()
###################
plt.subplot(2,3,4)
plt.title('Watts-Strogatz p = 0.1', fontsize=15)
k = int(av_degree/2)
p = 0.1
G = nx.watts_strogatz_graph(N, k, p, seed=None)
Tmax = 40
vlabda, rhoI_SIS, rhoI_SIR = critical_SIS_SIR_models(G, mu, Tmax)
plt.plot(vlabda, rhoI_SIS, 'ro--', label='SIS')
plt.plot(vlabda, rhoI_SIR, 'bo--', label='SIR')
plt.ylabel('Infected nodes', fontsize=13)
lbdac = moment_degree(G, 1)/moment_degree(G,2)
plt.axvline(lbdac, color='g', label='$\lambda_c$ = {0:5.3f}'.format(lbdac))
plt.legend()
####################
plt.subplot(2,3,5)
plt.title('Watts-Strogatz p = 0.01', fontsize=15)
k = int(av_degree/2)
p = 0.01
G = nx.watts_strogatz_graph(N, k, p, seed=None)
vlabda, rhoI_SIS, rhoI_SIR = critical_SIS_SIR_models(G, mu, Tmax)
plt.plot(vlabda, rhoI_SIS, 'ro--', label='SIS')
plt.plot(vlabda, rhoI_SIR, 'bo--', label='SIR')
plt.xlabel('$\lambda$', fontsize=13)
lbdac = moment_degree(G, 1)/moment_degree(G,2)
plt.axvline(lbdac, color='g', label='$\lambda_c$ = {0:5.3f}'.format(lbdac))
plt.legend()
####################
plt.subplot(2,3,6)
plt.title('Watts-Strogatz p = 0.001', fontsize=15)
k = int(av_degree/2)
p = 0.001
G = nx.watts_strogatz_graph(N, k, p, seed=None)
vlabda, rhoI_SIS, rhoI_SIR = critical_SIS_SIR_models(G, mu, Tmax)
plt.plot(vlabda, rhoI_SIS, 'ro--', label='SIS')
plt.plot(vlabda, rhoI_SIR, 'bo--', label='SIR')
plt.ylabel('Infected nodes', fontsize=13)
lbdac = moment_degree(G, 1)/moment_degree(G,2)
plt.axvline(lbdac, color='g', label='$\lambda_c$ = {0:5.3f}'.format(lbdac))
plt.legend()
plt.show()
We can see the presence of a phase transition point at the critical ratio probability $\lambda_c$ that increases from BA, ER and WS model, respectively. The lowest value belonging to ER and BA (they are very close to each other) can be explained by its densely random connected nodes (for ER), which provides a network enough connected to keep the infection self-sustainable for lower ratio of $\lambda$ and, for BA, the fact that hubs can spread easily the infection and, after the central node at the hubs recovery itself (we are considering only SIS model), its highly connected characteristic made it once again easily infected.
On the contrary, despite of create shortcuts on network's paths (small world property), its low average shortest path leghts allows the infection to be spreaded out, however, due to the close-to-average nodes' degrees (distinctly from BA), it demands a highly infectious infection to self-sustainability, which increases the value of $\lambda_c$
Question 3¶
Finally, we can study how the spread of infections take place on real networks and how do the nodes' properties (such as degree, local clustering coefficient, random walk accessibility, etc.) from which we place the first infected agent alters the infectious dynamics
def build_transit_probab_matrix(G):
"""
FUNCTION:
Build the Transition Probability Matrix
PARAMETERS:
G: networkx graph object
RETURN:
numpy matrix 2x2: Transition probability matrix address
"""
#Firstly, we can build a null vector in order to store
#each element from our probability matrix
adj_matrix = nx.to_numpy_matrix(G, nodelist=list(G.nodes()))
total_nodes = len(adj_matrix)
P_matrix = np.zeros((total_nodes, total_nodes))
#These loop will set the principal diagonal as zero
#and divide all other elements byt the sum of each
#row (normalizing each of them)
for i in range(total_nodes):
for j in range(total_nodes):
if adj_matrix[i,j] == 0:
P_matrix[i,j] = 0
else:
P_matrix[i,j] = adj_matrix[i,j]/np.sum(adj_matrix[i])
return P_matrix
def accessibility(G):
"""
FUNCTION:
Calculates the nodes accessibility of a graph ia
matrix exponential
PARAMETERS:
G: networkx graph object
RETURN:
"""
import scipy.linalg as alg
N = len(G.nodes())
P = build_transit_probab_matrix(G)
P2 = alg.expm(P)/np.exp(1)
vacc = np.zeros(N, dtype=float)
for i in np.arange(0, N):
acc = 0
for j in np.arange(0, N):
if P2[i,j] > 0:
acc += P2[i,j]*np.log(P2[i,j])
acc = np.exp(-acc)
vacc[i] = acc
return vacc
def degree_distribution(G):
"""
FUNCTION:
Calculates the degree distribution and the probability associated with each data
PARAMETERS:
G: networkx graph object
RETURN:
tuple (kvalues, pk): kvalues: each line segmentation for histogram-plot purpose;
pk: probability associated with each line segmentation
"""
vk = dict(G.degree())
vk = list(vk.values())
maxk = np.max(vk)
mink = np.min(min)
kvalues = np.arange(0, maxk+1)
pk = np.zeros(maxk+1)
for k in vk:
pk[k] = pk[k] + 1
pk = pk/sum(pk)
return kvalues, pk
def SIR_relate_node_metric(G, n, beta, mu, Tmax):
"""
FUNCTION:
Performs the SIR dynamic process on a given graph, placing a infectious agent (seed)
on each node each time and spreading it over the others
PARAMETERS:
G: networkx graph object
beta: probability of infection of a Susceptible agent
mu: probability of recovery of a Infected agent
Tmax: Maximum time over which the simulation will proceed
RETURN:
tuple (vt, av_infect): Contains vectors vt: time span
av_infect: average number of infectious nodes
"""
N = len(G)
all_infected = np.zeros(N)
all_recovered = np.zeros(N)
for k in range(n):
count = 0
for seed_node in G.nodes():
vector_states = np.zeros(N)
vector_states[seed_node] = 1
ninfected = 1
t = 0
infect = 0
recov = 0
while ninfected > 0 and t < Tmax:
infected = find(vector_states, 1)
for i in infected:
neigs = list(G.neighbors(i))
#print(neigs)
for j in neigs:
if np.random.rand() < beta:
if (vector_states[j] != 2):
vector_states[j] = 1
for k in infected:
if np.random.rand() < mu:
vector_states[k] = 2
ninfected = len(find(vector_states, 1))
nrecover = len(find(vector_states, 2))
infect+=ninfected/N
recov+=nrecover/N
t+=1
all_infected[count] += infect/n
all_recovered[count] += recov/n
count+=1
K = dict(G.degree())
K = list(K.values())
BC = dict(nx.betweenness_centrality(G))
BC = list(BC.values())
CC = dict(nx.closeness_centrality(G))
CC = list(CC.values())
EC = dict(nx.eigenvector_centrality(G))
EC = list(EC.values())
AC = accessibility(G)
KC = dict(nx.core_number(G))
KC = list(KC.values())
COMC = dict(nx.subgraph_centrality(G))
COMC = list(COMC.values())
df = pd.DataFrame({'Recovered':all_recovered,'Infected':all_infected,'K':K,'BC':BC,'CC':CC,'EC':EC,'AC':AC,'KC':KC, 'COMC':COMC})
print(df.head())
return df
def plot_matrix_correl(dataFrame):
"""
FUNCTION:
Plot the correlation matrix for several quantities
PARAMETERS:
dataFrame: pandas-library dataFrame object
RETURN:
None
"""
corr = dataFrame.corr()
plt.figure(figsize=(8,8))
plt.imshow(corr, cmap='Blues', interpolation='none',
aspect='auto')
plt.colorbar()
plt.xticks(range(len(corr)), corr.columns, rotation='vertical',
fontsize=13)
plt.yticks(range(len(corr)), corr.columns, fontsize=13)
plt.suptitle('Correlation between Centrality Measures',
fontsize=15)
plt.grid(False)
plt.show(True)
def plot_corr_metric_recovered(dataFr, suptitle):
plt.figure(1, figsize=(16,15))
plt.suptitle(suptitle, fontsize=20)
plt.subplot(3,3,1)
plt.scatter(dataFr['K'], dataFr['Recovered'], c='black')
plt.xlabel('K')
plt.subplot(3,3,2)
plt.scatter(dataFr['BC'], dataFr['Recovered'], c='black')
plt.xlabel('BC')
plt.subplot(3,3,3)
plt.scatter(dataFr['CC'], dataFr['Recovered'], c='black')
plt.xlabel('CC')
plt.subplot(3,3,4)
plt.ylabel('Fraction of recovered', fontsize=13)
plt.scatter(dataFr['EC'], dataFr['Recovered'], c='black')
plt.xlabel('EC')
plt.subplot(3,3,5)
plt.scatter(dataFr['AC'], dataFr['Recovered'], c='black')
plt.xlabel('AC')
plt.subplot(3,3,6)
plt.scatter(dataFr['KC'], dataFr['Recovered'], c='black')
plt.xlabel('KC')
plt.subplot(3,3,7)
plt.scatter(dataFr['COMC'], dataFr['Recovered'], c='black')
plt.xlabel('COMC')
plt.legend()
plt.show()
pos_hamst, G_hamster = read_txt_and_Gcc('/home/jorge/Documents/redes_complexas/petster-friendships-hamster/out.petster-friendships-hamster-uniq')
pos_celegans7, G_celegans7 = read_txt_and_Gcc('/home/jorge/Documents/redes_complexas/wi2007.txt', comment='#', nodetype=str)
pos_usair, G_usair = read_txt_and_Gcc('/home/jorge/Documents/redes_complexas/US_largest500_airportnetwork.txt')
Tmax = 20
beta = 0.3
mu = 1
n = 30
dataFr = SIR_relate_node_metric(G_celegans7, n, beta, mu, Tmax)
plot_corr_metric_recovered(dataFr, 'C. elegans neural network')
Taking the average over $n = 30$ iterations, we can see that the random walk accessibility (AC) has the best correlation between the fraction of recovered individuals. That is, starting the infection from a node with high accessibility, the higher the number of recovered nodes that will exist on the network after a given time span. Once the existence of recovered nodes implies they were infected (since they do not die, i.e. are removed from the network), we can conclude that the higher the accessibility of the initial infected node, the wider the infection spreads on C. elegans neural network.
Tmax = 20
beta = 0.3
mu = 1
n = 30
dataFr = SIR_relate_node_metric(G_hamster, n, beta, mu, Tmax)
plot_corr_metric_recovered(dataFr, 'Hamster petshop network')
For hamsters petshop friends, however, the metric that most correlates with the infection spreading is the Clustering coefficient.
Tmax = 20
beta = 0.3
mu = 1
n = 30
dataFr = SIR_relate_node_metric(G_usair, n, beta, mu, Tmax)
plot_corr_metric_recovered(dataFr, 'US airport network')
Finally, for US airport network, the highest correlation can be draw both from Accessibility and Clustering Coefficient (even though CC seems to have a better correlation).
From all those above studied networks, we can conclude that AC and CC happens to be the most relevant characteristic of a node to predicts how well a infection will spread out over a network.
And for last considerations, let us see how assortativennes on BA model changes the infection spread over a network. Thus, using the Xulvi-Brunet-Sokolo algorithm already built and explained, we can starts from a highly dissortative BA network and ends up toa highly assortative one while spreading infections over them
plt.figure(1, figsize=(15,12))
plt.suptitle('SIS and SIR spread for Barabasi-Albert with Assortativity', fontsize=20)
N = 1000
av_degree = 10
m = int(av_degree/2)
G = nx.barabasi_albert_graph(N, m)
beta = 0.3
mu = 1
Tmax = 40
p = 1
T = 15000
count = 1
for i in range(T+1):
xulvi_brunet_sokolov_disassortative(G, p)
T = 10000
for i in range(T+1):
xulvi_brunet_sokolov_assortative(G, p)
assort = nx.degree_assortativity_coefficient(G)
if i == 0 or i == 500 or i == 5000 or i ==10000:
vt, infSIR = SIR_model_spread(G, beta, mu, Tmax)
vt, infSIS = SIS_model_spread(G, beta, mu, Tmax)
#print(count)
plt.subplot(2,2,count)
count+=1
fc = fcritical(G)
plt.title('Assortativity: {0:.4f}'.format(assort), fontsize=15)
plt.plot(vt, infSIR, 'bo-', label='SIR')
plt.plot(vt, infSIS, 'ro-', label='SIS')
plt.legend()
if i == 0 or i == 5000:
plt.ylabel('Fraction of infected', fontsize=13)
if i == 5000 or i == 10000:
plt.xlabel('Time', fontsize=13)
#print('Done!')
plt.show()
We can see that the fraction of infected slightly reduces continuously from dissortative to assortative (up-left to bottom-right figure) and, most insterestingly, it also diminishes the noise (in this case, noise is related to how much doesthe number of infected agents vary over time), which is probably due to the fact that highly dissortative network have a few nodes connected with a bunch of low-degree ones and, because of that, considering that at the time $t_1$ the central node is infected, at time $t_2$, several nodes will be infected. If we consider that the central node became susceptible again, it will most certainly infected again while part of the low-degree's changed to susceptible once again. This cicle between peripheral nodes infected-sucseptible and central sucseptible-infected makes the huge oscillation that can be seen on the first two up plots.
plt.figure(1, figsize=(15,12))
plt.suptitle('SIS and SIR spread for Watts-Strogatz (p=0.05) with Assortativity', fontsize=20)
N = 1000
av_degree = 10
k = int(av_degree/2)
p = 0.01
G = nx.watts_strogatz_graph(N, k, p, seed=None)
beta = 0.3
mu = 1
Tmax = 100
p = 1
T = 15000
count = 1
for i in range(T+1):
xulvi_brunet_sokolov_disassortative(G, p)
T = 10000
for i in range(T+1):
xulvi_brunet_sokolov_assortative(G, p)
assort = nx.degree_assortativity_coefficient(G)
if i == 0 or i == 500 or i == 5000 or i ==10000:
vt, infSIR = SIR_model_spread(G, beta, mu, Tmax)
vt, infSIS = SIS_model_spread(G, beta, mu, Tmax)
#print(count)
plt.subplot(2,2,count)
count+=1
fc = fcritical(G)
plt.title('Assortativity: {0:.4f}'.format(assort), fontsize=15)
plt.plot(vt, infSIR, 'bo-', label='SIR')
plt.plot(vt, infSIS, 'ro-', label='SIS')
plt.legend()
if i == 0 or i == 5000:
plt.ylabel('Fraction of Infected', fontsize=13)
if i == 5000 or i == 10000:
plt.xlabel('Time', fontsize=13)
#print('Done!')
plt.show()
